function plot_basis_points_2d(basis_type, flag)

dN = reference_basis_points_2d(basis_type);

if nargin == 1
    Nb = size(dN, 1);
else
    Nb = min(flag, size(dN,1));
end

switch basis_type
    case {"P1", "P1b", "P2"}
        gauss_points = [(1-sqrt(3/5))/2, (1+sqrt(3/5))*(1-sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))*(1+sqrt(3/5))/4;
            1/2, (1-sqrt(3/5))/4;
            1/2, 1/4;
            1/2, (1+sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1-sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1+sqrt(3/5))/4];
        gauss_weights = [100/324*(1+sqrt(3/5))/8, 5/81*(1+sqrt(3/5)), 100/324*(1+sqrt(3/5))/8, 5/81, 8/81, 5/81, 100/324*(1-sqrt(3/5))/8, 5/81*(1-sqrt(3/5)), 100/324*(1-sqrt(3/5))/8];
    case {"Q1", "Q1b", "Q2"}
        gauss_points = [-sqrt(3/5), -sqrt(3/5);
            -sqrt(3/5), 0;
            -sqrt(3/5), sqrt(3/5);
            0, -sqrt(3/5);
            0, 0;
            0, sqrt(3/5);
            sqrt(3/5), -sqrt(3/5);
            sqrt(3/5), 0;
            sqrt(3/5), sqrt(3/5)];
        gauss_weights =[25/81, 40/81, 25/81, 40/81, 64/81, 40/81, 25/81, 40/81, 25/81];
end

x = reshape(gauss_points(:,1), 3, 3);
y = reshape(gauss_points(:,2), 3, 3);

for i = 1:Nb
    f = reshape(dN{i}(:), 3, 3);
    surf(x, y, f);
    hold on
end

end